Critical level-statistics for weakly disordered graphene 
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In two dimensions chaotic level-statistics is expected for massless Dirac fermions in the presence 
of disorder. For weakly disordered graphene flakes with zigzag edges the obtained level-spacing 
distribution in the Dirac region is neither chaotic (Wigner) nor localized (Poisson) but similar to 
that at the critical point of the Anderson metal-insulator transition. The quantum transport in 
finite graphene can occur via critical edge states as in topological insulators, for strong disorder the 
Dirac region vanishes and graphene behaves as ordinary Anderson insulator. 
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alistically large system sizes to occur. In weakly disor- 
dered graphene flakes with zz edges for open boundary 
conditions (be) the energy levels for a wide disorder re- 
gion (W < W c in Fig.l) show unusual quantum chaos, 
they resemble a gas of solitons Iff, etc. The correspond- 
ing level-spacing distribution is neither Wigner (chaotic) 
nor Poisson (localized), is critical as at the 3D Ander- 
son transition (It|. The energy spectrum, obtained by 
numerical diagonalization under the constraint of open 
be, gives our main result: weakly disordered graphene 
flakes with zz edges in the Dirac region show a novel 
kind of quantum chaos (W < W c in Fig.l) with critical 



Graphene, a monolayer of graphite with carbons ar- 
ranged in a two-dimensional honeycomb lattice has at- 
tracted great interest since its pioneering fabrication [J- 
4]. Its high mobility, current-carrying capacity and 
thermal conduction makes it a good candidate to re- 
place silicon in future nanoelectronics[f|. Graphene is 
a zero gap semiconductor, its Fermi surface consists of 
two non-equivalent Dirac points with linear energy dis- 
persion in the Dirac region. The sublattice (chiral) 
symmetry of graphene guarantees the presence of the 
Dirac region where electrons behave as massless Dirac 
fermions, limits backscattering, introduces Andreev re- 
flection and Klein tunnelling, etcjgj. In order to con- 
trol its gapless nature graphene was folded with cyclic 
boundary conditions in the short axis to form long cylin- 
ders (nanotubesQ), or long narrow strips with open per- 
pendicular edges (nanoribbons[H Q). In nanoribbons if 
the lattice geometry contains zigzag (zz) edgesj8-ll| the 
sublattice symmetry, which is responsible for the Dirac 
equation, guarantees thepresence of E = edge states 
tied to the boundary || We have studied Anderson 
localizationflij in finite disordered graphene samples of 
various shapes known as graphene flakes [3], that is in 
confined lattices fabricated from graphene. 

In two-dimensional(2D) disordered systems, with time- 
reversal and spin-rotation, all states are localized even 
for very small disorder, while in 3D localization occurs 
only above the Anderson metal-insulator transition (l^. 
In 2D, however, the states have localization lengths which 
often exceed the system size. In other words, for weak 
disorder the states behave as diffusive, they extend from 
one end of the sample to the other, and the energy lev- 
els give chaotic (Wigner) level-spacing distribution (ljj. 
The observed level-repulsion and spectral rigidity are the 
main characteristics of quantum chaos, the Poisson ex- 
ponential law which denotes localization requires unre- 
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FIG. 1: The energy levels E vs. the strength of disorder W 
for a quarter circle graphene flake of 9691 sites. For W < W c 
the level-statistics is critical and for W > W c is localized 
(Poisson). The levels were unfolded by removing variations 
of the density of states p(E) averaged over 5000 realizations 
of disorder. In this case W c — 2.3. 



Our calculations are done for the nearest neighbour 
tight-binding Hamiltonian 
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Ci(cj) annihilates (creates) an electron at A or B site i of 
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the honeycomb lattice, the diagonal disorder e$ is con- 
stant in the range [-W/2, +W/2], W is its strength and 
< i, j > denotes nearest neighbours with hopping matrix 
elements jij — 1. For W — the band structure displays 
two non-equivalent valleys (act as pseudospin states) re- 
lated by time-reversal [j| @] and in the Dirac region (at 
long length scales) Eq.(l) can be replaced by the con- 
tinuous Dirac equation[6|. The short-range diagonal dis- 
order via ti causes intervalley scattering and mixes the 
two-valleys which leads to localization. The decoupling 
of valleys by breaking time-reversal symmetry occurs for 
a smooth long-range disordered potential, this involves 
scattering within a single valley and gives weak antilo- 
calization familiar from spin-orbit coupling [l8j|. The dis- 
order e, mixes valleys and destroys A,B sublattice sym- 
metry. 

We have cut the honeycomb lattice in various shapes, 
circular, stadium, square, etc. The shape of the flake 
turns out to be irrelevant for our study but the perimeter, 
as we shall see, it is not. In our computations the brick- 
wall lattice is mostly used, however we find no signifi- 
cant differences from the honeycomb lattice itself. The 
eigenvalues of Eq.(l) in the Dirac region are obtained via 
Lanczos numerical diagonalization, building a statistical 
ensemble of random H matrices for every W . For W = 
quantum chaos for stadium flakes and integrability for 
circular flakes is expected, while if the lattice geome- 
try contains zz edges the sublattice symmetry guaran- 
tees the presence of edge states. The edge states due to 
the symmetry of the bulk belong to one type (A or B) 
sublattice and they are protected against disorder by the 
lattice topology which is reflected in the Hilbert space 
structure. In the honeycomb lattice the edge states arise 
from quantum interference of an incoming wave from one 
bond which splits into two. In the presence of weak disor- 
der the edge states spread almost uniformly in the Dirac 
region[Ioj ■ For zero disorder from the tight-binding equa- 
tions of the bulk[l, Q their amplitude ^(m) is nonzero 
only for m = and they are localized at the boundary 
having a fractal dimension of 1. For nonzero W they pen- 
etrate into the bulk and their fractal dimension varies. 
These massless degrees of freedom at the edges for fi- 
nite W are remnants of the sublattice symmetry which 
is broken by the diagonal disorder of Eq.(l). 

In Fig. 2 the averaged density of states p(E) for var- 
ious W is shown. In the absence of disorder (W = 0) 
the zz edges contribute to zero energy, if they disappear 
the edge states also disappear. Their total number re- 
mains fixed for nonzero disorder [loj], e.g. for W = 1,2 
the majority is in the almost constant p(E) region of 
Fig. 2. Moreover, for weaker disorder W = 0.01,0.1 mini- 
bands and minigaps develop in the spectrum around the 
discrete W = states, this is due (inset of Fig. 2) to the 
small p(E) and the finite lattice. We have also computed 
N(E) which counts the number of states from E = + 
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FIG. 2: The averaged density of states p(E) for circular 
graphene flakes with disorder W obtained from more than one 
million eigenvalues of Eq.(l). For W — 1,2 the edge states lie 
close to E — where p(E) is nonzero and almost constant. 
For strong disorder (W = 5) their proportion becomes van- 
ishingly small while in the opposite limit of pure (W = 0) 
graphene: p(E) = ■^jg\E\ (continuous line) is linear. Inset: 
the averaged integrated density of states for weaker disorder 
W = 0.01,0.1 the N(E) = J ( f+p(E') dE' (the E = states 
are not included) shows minigaps due to low p(E) values, 
while for higher W the structure disappears. 



to E. In a log-log plot we find constant N(E)/E vs E 
which implies constant p(E) for very low E, for slightly 
higher E the p(E) decreases for W < W c (due to a mini- 
gap) and it increases for W > W c (no minigaps). For 
strong disorder (W = 5) states fill up the Dirac region 
and make p(E) constant. 

We have examined in detail the perimeter of the con- 
sidered flakes, identified each type of lattice edge and its 
contribution to p{E). We find armchair and zz edges (the 
dangling bonds arc rare), the most frequent edges arc the 
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The ratio of zz to armchair edges for circular 
flakes is ~ 3.8 and it varies linearly with the averaged 
flake radius R (the linear curve showed more oscillations 
for the zz edges). The ratio of zz edges over the total 
number of sites tends to zero inversely proportional to 
R, since their number is ~ R and the total number of 
sites is proportional to the area of the flake ~ R 2 , the 
density of states p(0) reaches a maximum before it van- 
ishes as 1/R for large R. The disorder shifts the edge 
states away from zero energy without changing their to- 
tal number (Io|. 

Let us now discuss the eigenvalues of Eq.(l). We ex- 
amined the level-statistics of the unfolded energy spec- 
tra via the nearest spacing distribution P(S) which 
can distinguish between chaotic (Wigner) and localized 
(Poisson) [14| . In Fig. 3 the obtained P(S) of the first two 
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infinite size and/or strong disorder unlike in topological 
in □ 



FIG. 3: The calculated spacing distribution P(S) of the first 
two positive energies (S = E2 — E\) for W < W c for circular 
quarter graphene flakes in the new phase W < W c , sizes N — 
2052, 9691, 15638 and 50000 realisations of disorder. The data 
are independent of size and are described by the intermediate 
semi-Poisson P(S) = 4S exp(—2S) distribution (black line). 



positive levels (S = E2 — E%) is shown for the critical 
region of W < W c . The P(S) is independent of size and 
fits into a curve which has mixed chaotic and localized 
character, it is chaotic for small spacing S and localized 
for large S"|l7], it is described by P(S) = 4S , exp(-2S*) 
and belongs to the semi-Poisson family [l9j]. The critical 
P(S) interpolates between chaotic and localized limits: 
P(S) ~ S for 5 ^ and P(S) ~ exp(-(l + /3)S) for 
S — > 00, the universality class index f3 = 1, for broken 
time-reversal symmetry via a magnetic field /3 = 2[l3| 
and for broken spin rotation via spin-orbit coupling /3 = 
4R1. 



For pure graphene (W — 0) we find Poisson level- 
statistics in the Dirac region since the integrated den- 
sity of states N(E ~ ^fk§ + k®) ~ E 2 gives integrable 

unfolded levels an 2 + m 2 , n,m — 1,2,..., a ~ 1 for a 
square sample. For graphene with zz edges the critical 
distribution replaces Poisson for W < VFc(Fig.l) which 
can describe only larger S, while for small S the P(S) is 
Wigner-like. The obtained critical distribution interpo- 
lates between chaotic and the localized limits (apart from 
3D critical sy_stems it was found for pseudo- integrable bil- 
liards, etc.[19J). In graphene localization is easier than 
other 2D due to its small coordination, this counterbal- 
ances the chiral direction of electrons in the Dirac re- 
gion which favours absence of localization. The effect 
of edges vanishes by taking periodic boundary condi- 
tions, as it was seen in previous numerical studies of 
level-statistics 2l|, |22[ . For strong disorder the Dirac re- 
gion vanishes altogether and Poisson is obtained for lo- 
calized states as for zero disorder. The flow of P(S) to- 
wards Poisson as the size increases also occurs for the 
honeycomb lattice but for a smaller W than for the 
square. The topological effects found in Fig.l vanish for 
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FIG. 4: The flow of the level-spacing distribution P(S) to- 
wards Poisson as the system size increases. The data for 
stronger disorder W > W c , W — 1,5 are obtained from 
stadium quarter flakes of sizes N — 924, 2967, and 8940, 
for 50000 realisations of disorder and energies in the window 
[0, 0.12] of the Dirac region. The two arrows indicate the flow 
towards the localized (Poisson) limit as the system size in- 
creases, the approach although slow is much faster than that 
for the square lattice. The Poisson limit verifies that in infi- 
nite graphene all states are localized. 



FIG. 5: A critical edge state at nonzero E = 0.14 for a quarter 
circular flake with disorder W = 2 (W < W c )- Its fractal 
dimension is close to 1. 



Our work was partially motivated by the fabrication 
of graphene quantum dots for diameters ranging from 
40 to 100 nanometers [i~3j . For small flakes the experi- 
mental level-statistics is increasingly non-Poisson while 
the very small flakes of a few nm width remained always 
conductive. The obtained critical level-statistics is differ- 
ent from that of conventional 2D disordered systems. We 
find that the sample size is also crucial, the fast or slow 
approach to localization depends on the kind of lattice, 
in graphene for strong disorder W > W c the approach 
to localization is faster than in other 2D (Fig. 4). For 
weak disorder W < W c the level-statistics depends on 
the lattice termination, e.g for a nanoribbon of arbitrary 
orientation which has a majority of zz edges [11] critical 
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level-statistics should also appear. The edge states are 
prominent only for small size and weak disorder and are 
intimately connected to the small p(E) in the Dirac re- 
gion. It is also interesting to ask whether the critical 
level-statistics can give conductive transport. 

The edge states in flakes and nanoribbons are pro- 
tected from disorder, they cannot be deleted without 
changing the topology of the Hilbert space. Edge states 
of topological origin are found for bilayer graphene with 
subgap conductance [24|. are seen in microwave real- 
izations of nanoribbons |25|, they are spatially resolved 
with scanning tunnelling microscopy [26j, etc. The new 
phase found for W < W c also shows the absence of 
weak localization[ll-0] in graphene. The boundary ge- 
ometry plays important role, the critical P(S) is inde- 
pendent of size but it depends on be as at the 3D An- 
derson transition (2?]]. In the presence of a magnetic field 
transport via protected edge states occurs in the integral 
quantum Hall effect fljj ■ It would be interesting to exam- 
ine our findings in the presence of a magnetic field, for 
other types of disorder, such as edge disorder, resonant 
impurities, vacancies, adding the true spin of electrons 
and electron-electron interactions. 

The level-statistics of weakly disordered graphene 
flakes as a function of short-ranged disorder is not chaotic 
as in other 2D where the localization length is usually 
much larger than the system size. The realistic edge 
structure for W < W c makes the spectrum in the Dirac 
region intermediate between chaotic and integrable, sim- 
ilar to that at the critical point of the Anderson metal- 
insulator transition. Graphene is different than other 
weakly disordered 2D, its topological behaviour is also 
suggested by the finite conductivity a ~ e 2 /h at the 
Dirac point [6j. On the basis of our results a mapping 
of electrons to the classically integrable Calogero-Moser 
model of interacting particles is possible [lij], a kink mech- 
anism relates chaotic behavior to classical equations. In 
graphene the edge is important and when combined with 
the low p(E) gives the novel critical regime of Fig.l. 
The smaller coordination number of the lattice affects 
a graphene flake, we find smaller and smaller W c by in- 
creasing the flake's size and in the infinite size limit the 
spectrum becomes gapless. In other words, graphene has 
no intrinsic topological order but a pseudogapped phase 
for weak disorder and finite size. 

In conclusion, graphene for weak disorder can conduct 
via its edges. A sharp transition at W c is found (Fig.l) 
which distinguishes between critical (W < W c ) and lo- 
calized (W > W c ) states. This is shown in the Dirac 
region where the density of states is low and minibands 
appear with edge states. Our results are a signature that 
weakly disordered graphene is topological insulator-like, 
its conduction depends on the topology of the perimeter. 
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